home *** CD-ROM | disk | FTP | other *** search
/ TPUG - Toronto PET Users Group / TPUG Users Group CD / TPUG Users Group CD.iso / PET / S-Super PET / (s)t2.d64 / PERIODIC.FTN < prev    next >
Encoding:
Text File  |  2009-01-18  |  12.6 KB  |  442 lines
  1. *              Program  periodicities
  2. *
  3. *  This program was typed directly from the book:
  4. *
  5. *    FOURIER ANALYSIS OF TIME SERIES: AN INTRODUCTION
  6. *    by Peter Bloomfield
  7. *    John Wiley & Sons, 1976
  8. *
  9. *  This main program and subroutine  datin  are used to
  10. *  input data for, and output results from, fitting a
  11. *  model of hidden periodicities to a time series.
  12. *  The time series data is read by  datin .  Data input
  13. *  from a control file are  conv , lim , and appflg  (which
  14. *  are defined in subroutine  optom ),  the number of
  15. *  frequencies to be fitted, followed by the list of
  16. *  approximate frequencies and the estimated sine and
  17. *  cosine coefficients.
  18. *
  19. *  Note:- An approximate fit is obtained if  appflg and lim
  20. *         are set to unity and conv is set to pi or similar
  21. *         large number.  An 'exact' fit is obtained with
  22. *         appflg set to zero, lim to 5 (for example) and
  23. *         conv set to a small number like 1.0e-05 to allow
  24. *         the algorithm to iterate to convergence.
  25. *
  26. *
  27.       real x(600),y(600),fre(10),a(10),b(10)
  28.       integer appflg
  29.       character fname1,fname2,fname3
  30. *
  31.       write(6,1)
  32.     1 format("1Enter the time series file name")
  33.       read,fname1
  34.       print,"Enter the control parameter filename"
  35.       read,fname2
  36.       print,"Enter the output file name"
  37.       read,fname3
  38. *
  39.       call datin (x,n,start,step,1,fname1)
  40.       open(unit=2,file=fname2)
  41.       read(2) conv,lim,appflg
  42.       open(unit=3,file=fname3)
  43.       write(3,4) conv,lim,appflg
  44.     4 format(" For this run,  conv =",e12.4/
  45. &            "                 lim =",i5/
  46. &            "              appflg =",i5)
  47.       read(2) m
  48.       read(2) (fre(j),a(j),b(j),j=1,m)
  49.       write(3,5)
  50.     5 format(/" Initial values are -")
  51.       write(3,3) (j,fre(j),a(j),b(j),j=1,m)
  52.     3 format(" Component Frequency     Cosine         Sine"/
  53. &            "                             Coefficients"/
  54. &            (i5,f15.7,2e15.6))
  55. *
  56.       call optom (x,n,rmu,fre,a,b,m,conv,lim,appflg,y)
  57. *
  58.       write(3,6)
  59.     6 format(/"  Final values are -")
  60.       write(3,3) (im,fre(im),a(im),b(im),im=1,m)
  61.       ss = 0.0
  62.       do 30 i = 1 , n
  63.          temp = x(i) - rmu
  64.          do 40 j = 1 , m
  65.             arg = float(i-1) * fre(j)
  66.             temp = temp - a(j) * cos(arg) - b(j) * sin(arg)
  67.    40    continue
  68.          ss = ss + temp**2
  69.    30 continue
  70.       write(3,7) rmu,ss
  71.     7 format(/" Fitted constant is        ",e15.6/
  72. &             " Residual Sum of Squares is",e15.6)
  73.       stop
  74.       end
  75.  
  76.  
  77.       subroutine datin ( x , n , start , step , m , fname)
  78. *
  79. *  This subroutine is used to input a series of values
  80. *  (in run time format) and some associated quantities
  81. *  (in fixed format).  The first four data records are -
  82. *
  83. *  1   A header record (72 bytes)
  84. *  2   Value of n (i5)
  85. *  3   The data format (72 bytes)
  86. *  4   Start and step  (2f10.5)
  87. *
  88. *  Parameters are -
  89. *
  90. *  Name      Type                     Value
  91. *                         On Entry                On Return
  92. *
  93. *  x      Real array      Not used              The series
  94. *
  95. *  n      Integer         Not used              Series length
  96. *
  97. *  start  Real            Not used              Time value at the
  98. *                                               first data point
  99. *
  100. *  step   Real            Not used              Time increment
  101. *                                               between data points
  102. *
  103. *  m      Integer         Logical unit number   Unchanged
  104. *
  105. *  fname  Character       input file name      Unchanged
  106. *
  107.       real x(512)
  108.       character head , form , fname
  109.       open(unit=m,file=fname,access="sequential",recl=80)
  110.       read(m,1) head,n,form,start,step
  111.     1 format(a72/i5/a72/2f10.5)
  112.       write(6,2) head,n,form,start,step
  113.     2 format("0The data header reads -"/1x,a72/
  114. &            " The series length is",i6/
  115. &            " The data format is -"/1x,a72/
  116. &            " time origin is",f11.5,
  117. &            ",  time increment is",f11.5)
  118.       read(m,form) (x(i),i=1,n)
  119.       close(unit=m)
  120.       return
  121.       end
  122.  
  123.  
  124.       subroutine optom (x,n,rmu,fre,a,b,m,conv,lim,appflg,y)
  125. *
  126. *  This subroutine, with subprograms  localm ,  ssreg ,
  127. *  stats  and  parms , implements the algorithm for
  128. *  the least-squares fitting of the model of hidden
  129. *  periodicities.  Parameters are:-
  130. *
  131. *  Name     Type                         Value
  132. *                           On Entry                On Return
  133. *
  134. *  x      Real array    The time series           Unchanged
  135. *
  136. *  n      Integer       Series length             Unchanged
  137. *
  138. *  rmu    Real          Not used                  Constant term
  139. *
  140. *  fre    Real array    Starting values for       Final values
  141. *                       the frequencies to
  142. *                       be fitted
  143. *
  144. *  a      Real array    Starting values for       Final values
  145. *                       cosine coefficients
  146. *
  147. *  b      Real array    Starting values for       Final values
  148. *                       sine coefficients
  149. *
  150. *  m      Integer       Number of components      Unchangd
  151. *                       to be fitted
  152. *
  153. *  conv   Real          Convergence criterion     Unchanged
  154. *                       Iteration ceases when
  155. *                       in one cycle,  no
  156. *                       frequency changes by
  157. *                       more than conv
  158. *
  159. *  lim    Integer       Maximum number of         Unchanged
  160. *                       cycles of iteration
  161. *
  162. *  appflg  Integer      Flag controlling          Unchanged
  163. *                       whether approximate
  164. *                       ( 1 ) or exact ( 0 )
  165. *                       least squares is to
  166. *                       be usd
  167. *
  168. *  Notes.
  169. *  1.  If values of  n  or  m  exceeding 600 and 10.
  170. *      respectively, are used, the dimension statement
  171. *      below should be changed accordingly.
  172. *  2.  Use of exact least squares (appflg = 0) will
  173. *      cause longer execution times.
  174. *  3.  The frequencies in general converge to values
  175. *      within  2*pi/n  of their starting values.  The
  176. *      starting values should be given to at least this
  177. *      accuracy, or the algorithm may be trapped by
  178. *      side-lobes.  Starting values of the cosine and
  179. *      sine coefficients are less critical, and may be
  180. *      set to zero.
  181. *
  182.       real localm
  183.       real x(n) , y(n) , fre(10) , a(10) , b(10)
  184.       integer appflg
  185. *
  186. *     *** root of relative machine epsilon
  187. *
  188.       eps = 1.0
  189.     1 eps = 0.5 * eps
  190.       epsp1 = eps + 1.0
  191.       if (epsp1 .gt. 1.0) go to 1
  192.       eps = sqrt(eps)
  193. *
  194.       pi = 4.0 * atan(1.0)
  195.       t  = conv
  196.       delta = pi / float(n)
  197.       do 10 kount = 1 , lim
  198.          sum = 0.0
  199.          do 20 i = 1 , n
  200.             y(i) = x(i)
  201.             do 30 j = 1 , m
  202.                arg = float(i-1)* fre(j)
  203.    30       y(i) = y(i) - a(j) * cos(arg) - b(j) * sin(arg)
  204.    20    sum = sum + y(i)
  205.          rmu = sum / float(n)
  206.          test = 0.0
  207.          do 40 j = 1 , m
  208.             do 60 i = 1 , n
  209.                y(i) = x(i) - rmu
  210.                do 50 k = 1 , m
  211.                   if (k .eq. j) go to 50
  212.                   arg = float(i-1) * fre(k)
  213.                   y(i) = y(i) -a(k) * cos(arg) - b(k) * sin(arg)
  214.    50          continue
  215.    60       continue
  216.             dummy = localm (fre(j) - delta , fre(j) + delta,
  217. &                           eps , t , temp , y , n , appflg)
  218.             test = amax1 (test , abs(fre(j)-temp))
  219.             fre(j) = temp
  220.    40    call parms (y , n , fre(j) , appflg , a(j) , b(j))
  221.          if (test .lt. conv) return
  222.          delta = test + 2.0 * t
  223.    10 continue
  224.       return
  225.       end
  226.  
  227.  
  228.       real function localm (a , b , eps , t , x , y , n , appflg)
  229. *
  230. *  This is the FORTRAN function LOCALM given by
  231. *  Richard Brent in
  232. *       Algorithms for minimization without derivatives
  233. *  (Prentice-Hall, 1973).
  234. *  It finds a local minimum of the function  f  in the
  235. *  interfal (a , b).
  236. *  t  and  eps  define a tolerance  tol = eps * abs(x) + t ,
  237. *  where  x  is the current approximation to the position
  238. *  of the minimum.  The minimum is found with an error
  239. *  of at most  3 * tol.
  240. *  f  is not evaluated at points closer than  tol.
  241. *  a suitable value for  eps  is the square root of the
  242. *  relative machine precision.  For more detail see the
  243. *  above reference.
  244. *
  245.       real m , y(n)
  246.       integer appflg
  247.       sa = a
  248.       sb = b
  249.       x = sa + 0.381966 * (sb - sa)
  250.       w = x
  251.       v = w
  252.       e = 0.0
  253. *     fx = f(x)
  254.       fx = -ssreg (y , n , x , appflg)
  255.       fw = fx
  256.       fv = fw
  257.    10 m = 0.5 * (sa + sb)
  258.       tol = eps * abs(x) + t
  259.       t2 = 2.0 * tol
  260.       if (abs(x-m) .le. t2-0.5*(sb-sa)) go to 190
  261.       r = 0.0
  262.       q = r
  263.       p = q
  264.       if (abs(e) .le. tol) go to 40
  265.       r = (x - w) * (fx - fv)
  266.       q = (x - v) * (fx - fw)
  267.       p = (x - v) * q - (x - w) * r
  268.       q = 2.0 * (q - r)
  269.       if (q .le. 0.0) go to 20
  270.       p = -p
  271.       go to 30
  272.    20 q = -q
  273.    30 r = e
  274.       e = d
  275.    40 if (abs(p) .ge. abs(0.5*q*r)) go to 60
  276.       if ((p .le. q*(sa-x)) .or. (p .ge. q*(sb-x))) go to 60
  277.       d = p / q
  278.       u = x + d
  279.       if ((u-sa .ge. t2) .and. (sb-u .ge. t2)) go to 90
  280.       if (x .ge. m) go to 50
  281.       d = tol
  282.       go to 90
  283.    50 d = -tol
  284.       go to 90
  285.    60 if (x .ge. m) go to 70
  286.       e = sb - x
  287.       go to 80
  288.    70 e = sa - x
  289.    80 d = 0.381966 * e
  290.    90 if (abs(d) .lt. tol) go to 100
  291.       u = x + d
  292.       go to 120
  293.   100 if (d .le. 0.0) go to 110
  294.       u = x + tol
  295.       go to 120
  296.   110 u = x - tol
  297. * 120 fu = f(u)
  298.   120 fu = -ssreg(y , n , u , appflg)
  299.       if (fu .gt. fx) go to 150
  300.       if (u .ge. x) go to 130
  301.       sb = x
  302.       go to 140
  303.   130 sa = x
  304.   140 v = w
  305.       fv = fw
  306.       w = x
  307.       fw = fx
  308.       x = u
  309.       fx = fu
  310.       go to 10
  311.   150 if (u .ge. x) go to 160
  312.       sa = u
  313.       go to 170
  314.   160 sb = u
  315.   170 if ((fu .gt. fw) .and. (w .ne. x)) go to 180
  316.       v = w
  317.       fv = fw
  318.       w = u
  319.       fw = fu
  320.       go to 10
  321.   180 if ((fu .gt. fv) .and. (v .ne. x) .and. (v .ne. w)) go to 10
  322.       v = u
  323.       fv = fu
  324.       go to 10
  325.   190 localm = fx
  326.       return
  327.       end
  328.  
  329.  
  330.       function ssreg (y , n , omega , appflg)
  331. *
  332. *  This function returns the sum of squares (approximate
  333. *  or exact) associated with omega.  Parameters are
  334. *
  335. *  Name      Type          Value on Entry (none are changed)
  336. *
  337. *  y        Real array    The time series
  338. *
  339. *  n        Integer       Series length
  340. *
  341. *  omega    Real          The frequency
  342. *
  343. *  appflg   Integer       1   for approximate least squares
  344. *                         0   for exact least squares
  345. *
  346.       real y(n)
  347.       integer appflg
  348.       call stats (y , n , omega , cy , sy)
  349.       if (appflg .eq. 1) go to 10
  350.       rn = n
  351.       con = sin(rn*omega)/sin(omega)
  352.       arg = (rn-1.0)*omega
  353.       cc = 0.5*(rn+cos(arg)*con)
  354.       cs = 0.5*sin(arg)*con
  355.       ss = rn-cc
  356.       ssreg = (ss*cy**2-2.0*cs*cy*sy+cc*sy**2)/(cc*ss-cs**2)
  357.       return
  358.    10 continue
  359.       ssreg = (cy**2+sy**2)*2.0/float(n)
  360.       return
  361.       end
  362.  
  363.  
  364.       subroutine stats (y , n , omega , cy , sy)
  365. *
  366. *  This subroutine returns the cosine and sine sums of a
  367. *  time series.  Parameters are
  368. *
  369. *  Name      Type                       Value
  370. *                             On Entry         On Return
  371. *
  372. *  y       Real array    The time series     Unchanged
  373. *
  374. *  n       Integer       Series length       Unchanged
  375. *
  376. *  omega   Real          The frequency       Unchanged
  377. *
  378. *  cy      Real          Not used            Cosine sum
  379. *
  380. *  sy      Real          Not used            Sine sum
  381. *
  382.       real y(n)
  383.       cy = 0.0
  384.       sy = 0.0
  385.       do 10 i = 1 , n
  386.          arg = float(i-1) * omega
  387.          cy  = cy + cos(arg) * y(i)
  388.          sy  = sy + sin(arg) * y(i)
  389.    10 continue
  390.       return
  391.       end
  392.  
  393.  
  394.       subroutine parms (y , n , omega , appflg , a , b)
  395. *
  396. *  This subroutine returns the (exact or approximate)
  397. *  least squares estimate of the cosine and sine
  398. *  coefficients of a single periodic component.
  399. *  paramters are
  400. *
  401. *  Name       Type                       Value
  402. *                             On Entry           On Return
  403. *
  404. *  y        Real array   The time series        Unchanged
  405. *
  406. *  n        Integer      Series length          Unchanged
  407. *
  408. *  omega    Real         The frequency          Unchanged
  409. *
  410. *  appflg   Integer      1  or  0  for          Unchanged
  411. *                        approximate or
  412. *                        exact least squares,
  413. *                        respecively
  414. *
  415. *  a        Real         Not used               cos coefficient
  416. *
  417. *  b        Real         Not used               sin coefficient
  418. *
  419.       real y(n)
  420.       integer appflg
  421.       call stats (y , n , omega , cy , sy)
  422.       rn = float(n)
  423.       if (sin(omega) .eq. 0.0) go to 20
  424.       if (appflg .eq. 1) go to 10
  425.       con = sin(rn*omega)/sin(omega)
  426.       arg = (rn-1.0)*omega
  427.       cc = 0.5*(rn+cos(arg)*con)
  428.       cs = 0.5*sin(arg)*con
  429.       ss = rn-cc
  430.       del = cc*ss-cs**2
  431.       a = (cy*ss-sy*cs)/del
  432.       b = (sy*cc-cy*cs)/del
  433.       return
  434.    10 continue
  435.       a = 2.0*cy/rn
  436.       b = 2.0*sy/rn
  437.       return
  438.    20 a = cy/rn
  439.       b = 0.0
  440.       return
  441.       end
  442.